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Abstract 



We present a self-consistent electronic structure calculation method 
based on the Exact Muffin-Tin Orbitals (EMTO) Theory developed by 

0. K. Andersen, O. Jepsen and G. Krier (in Lectures on Methods of 
Electronic Structure Calculations, Ed. by V. Kumar, O.K. Andersen, 
A. Mookerjee, Word Scientific, 1994 pp. 63-124) and O. K. Andersen, 
C. Arcangeli, R. W. Tank, T. Saha-Dasgupta, G. Krier, O. Jepsen, and 

1. Dasgupta, (in Mat. Res. Soc. Symp. Proc. 491, 1998 pp. 3-34). 
The EMTO Theory can be considered as an improved screened KKR 
(Korringa-Kohn-Rostoker) method which is able to treat large overlap- 
ping potential spheres. Within the present implementation of the EMTO 
Theory the one electron equations are solved exactly using the Green's 
function formalism, and the Poisson's equation is solved within the Spher- 
ical Cell Approximation (SCA). To demonstrate the accuracy of the SCA- 
EMTO method test calculations have been carried out. 
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I. INTRODUCTION 



During the last decades many attempts have been made to develop accurate and at 
the same time efficient methods for solving the Kohn-Sham equations in an application 
of the Density Functional Theory for condensed matter. The accuracy of the methods 
is crucial e.g. when one searches for the answers given by different density functional 
approximations. The full-potential techniques have been specially designed to fulfill this 
requirement. Though, in principle, they give highly accurate results, they have their 
own limitations. In many cases a compromise has been made between the accuracy 
and efficiency, and methods based on approximate one electron potentials have been 
developed. The most commonly used muffin tin approach, albeit its mathematical for- 
mulation is very elegant, presents a rather poor representation of the exact potential. 
Though, the Atomic Sphere Approximation (ASA) brings a real improvement to the 
potential, most of the conventional methods based on the ASA use similar approxima- 
tion to the one electron energies and charge density as well ]I]. Therefore, using these 
methods, reasonably accurate results can only be obtained for close packed systems, 
and they are not suitable to treat systems of low symmetry. In order to maintain or 
increase the accuracy different corrections should be included and, therefore, the ASA 
based methods lose their elegance and efficiency. 

A few years ago breakthrough was made by developing the Exact Muffin-Tin Or- 
bitals (EMTO) Theory Within the EMTO Theory the one electron states are 

calculated exactly for the overlapping muffin-tin potential, while the solution of Pois- 
son's equation can include certain shape approximations, if required. By separating 
the two approaches used for the potential and one electron states the accuracy can be 
sustained at a level comparable with that of the full-potential techniques without losing 
significantly from the efficiency. The EMTO Theory can be considered as an improved 
screened KKR method [|],[5| , within that large overlapping potential spheres can be used 
for accurate representation of the exact one electron potential |3]||. 

In this work we present a self-consistent implementation of the EMTO Theory within 
the Spherical Cell Approximation (SCA) for the Poisson's equation. In the first part we 
review the EMTO Theory [@||, the definition of the screened spherical waves and the 
matching equation. Furthermore, we establish the expressions for the number of states 
and electron density using the Green's function formalism. In the second part of the 
paper we discuss the impact of the SCA, used for the shape of the Wigner-Seitz cell, 
on the total charge density and on the overlapping muffin tin potential. Finally, we 
establish the accuracy of the SCA-EMTO method by performing test calculations for 
some systems where reliable full-potential data are available. An approximate solution 
of the kink cancellation equation in order to reduce the number of iterations needed in 
a self-consistent calculation is presented in the Appendix. 

II. OVERVIEW: THE EMTO THEORY 

In the following we review the basic concepts of the EMTO Theory developed by 
O. K. Andersen and co-workers [0,0,0- Assume that the one-electron Kohn-Sham 
equations, 
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-V 2 + v(r)\ ^(r) = t 3 ^(r), (1) 
are solved within the muffin tin approximation for the effective potential, 

v(r) « v + KM - v ] , (2) 

where R runs over the lattice sites. Here and in the following we use the notation 
tr = r — R and omit the vector notation for the index R. In (Q) v o denotes a parameter 
that reduces to the muffin tin zero in the case of non-overlapping muffin tins. The 
spherical potentials Vr(tr) become equal to vq outside the potential spheres of radii 
sr. These radii can be chosen as the linear overlap between the spheres to be as large 
as 30 — 40 % [H|,|3]||. It has turned out that for a good representation of the real 
full-potential in terms of overlapping muffin tin wells usually a big overlap is preferred 
between the potential spheres ||. 

In order to solve the Schrodinger equation ([I]) for the muffin tin potential (0) one 
chooses different basis functions inside the potential spheres and in the interstitial 
region. Inside the sphere at R the partial waves are chosen as the basis function which 
are defined as the products of the regular solutions of the radial Schrodinger equation, 



d 2 



tr 0ra(e,r R ) 



1(1 + 1) 



r 2 R 



+ v R (r R )-e r R (j> m {e,r R ), (3) 



dr R 

and the real spherical harmonics, viz. 

(Prl^^r) = <l)Ri(e,r R )Y L (f R ), (4) 

where L = (I, m). Outside the spheres the so called screened spherical waves, ip RL (n, r R ), 
are used as basis functions. Therefore, the wave function for the energy e 3 - can be written 

as 

= E Me jt T R ) e R (r R ) u a RL>j + £ ^ a RL (^ h r R ) [1 - Q r (t r )] v a RLJ . (5) 

RL RL 

Here k 2 = e — vq, and in the non-overlapping muffin tins limit it denotes the interstitial 
one-electron kinetic energy. The 0_r(i"/j) is one inside the sphere of radius s R centered 
at R and zero outside. The expansion coefficients u RL j and v RL j as well as the en- 
ergies 6j are determined from the condition that the wave function ^(r) and its first 
derivative should be continuous at the potential spheres. The algebraic formulation of 
this matching condition in the EMTO formalism is the so the called kink cancellation 
equation, which is equivalent to the KKR (Korringa-Kohn-Rostoker) equation in an 
arbitrarily screened representation ||. 



A. The screened spherical waves 

The screened spherical waves can be defined in conjunction with hard spheres 
centered at all sites R with radii a^. They are solutions of the wave equation, 
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V 2 +k 2 r RL (K 2 ,r R ) 



0. 



(6) 



with the boundary condition that on their own a— spheres they behave like a pure real 
spherical harmonic, while the Yy{r R i) projections on all the other a— spheres, R' 7^ R, 
vanish. They form a complete basis set in the "a" interstitial region and may be 
expressed in terms of the "value", f RL , and the "slope", g RL , functions @, whose radial 
part satisfy the following boundary conditions 



r )\a r 



= 1 and 
and 



dfm( K > r ) 
dr 



Or 



\am 0' 
1 

11 ~ ~ ' 



(7) 
(8) 



These functions may of course be expressed in terms of the usual spherical Bessel and 
Neumann functions [|J 



fm(^r)=ji(Kr)W a {f,Kn} - k rn(Kr)W a {f 3} 



and 



9mi^ r ) = 3i(Kr)W a {g,Kn} - Kni(Kr)W a {g,j} 
since W r {j, n} — 1/k, and therefore satisfy the Wronskian 



dr 



dr 



am- 



(9) 



(10) 



The screened spherical wave V^zX^ 2 ? r n) rnay be expanded in spherical harmonics 
Yl'{tr') about any site R', as 



^rl(^^r) = fm(^r R )Y L (f R )5 RRI 5 W + Y.9m'{^ r w) Y L >(r R> ) S r , lirl (k 2 

L> 



:i2) 



where the expansion coefficients, S r , l , rl (k 2 ), are the elements of the so called slope 
matrix. The slope matrix can be derived from the bare KKR structure constant matrix 
B r >l',rl(k), by matrix inversion p| 



S a (K 2 ' 



V{j{K,a)} 



1 



-B(k) + k cot oi(k) 



-1 1 



j(K,a) 



(13) 



where T> denotes the logarithmic derivative, T>{j(r)} = r [dj(r)/dr] /j(r), and for sim- 
plicity where we have used matrix notation. We note that this equation is equivalent 
to Eq. (3.26) from Ref. and Eq. (15) from Ref. B. The bare KKR structure con- 
stants are defined as the expansion coefficients of the k til(k, r R ) = k rii(Kr R ) Yl{t r ) 
functions around site R' in terms of the jz/( K > r R') = 3i'( Kr R') Yyif^) functions, i.e. 

Kn L (K,r R ) = J2jL>(K,r R ,) B r , l , rl (k), 



with 
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B rilirl {k) = ^Y. c ll> rl+l '~ l " K ni/ ,(K,R'-R), (14) 

L" 

and where C^' L , are the real Gaunt numbers. 

For the partial waves explicitly included in the formalism, the so called low partial 
waves with I < l\ ow = 2 — 3, «#/(«;) are the hard sphere phase shifts given by 

cotaju(K) = rii(Ka m )/ji(Ka m ) 

and for the remaining i?/-chanels, a R i(n) are the proper phase shifts. For high /'s the 
latter vanish, and at that point the matrix to be inverted in (|13D can be truncated. 

When the hard sphere radii, a R , are properly chosen and k 2 lies below the bottom of 
the hard sphere continuum, the screened spherical waves have short range. Therefore, 
the slope matrix can be calculated in real space and the method is suitable to treat 
impurities, defects, surfaces, etc. It was shown in Ref. [0 that the shortest range of 
the screened spherical waves can be achieved for non-overlapping spheres with a R w 
0.5 — 0.85 s R , depending on the maximal orbital quantum number I of the partial waves 
explicitly included in the formalism. The s R denotes the inscribed or touching sphere 
radii. In the KKR community, it is customary to determine the a^^O's as the phase 
shifts of repulsive potentials. 

Because a screened spherical wave has pure (/, m) character only on its own a-sphere, 
the matching condition in Eq. (|5]) should be set up at this sphere. The connection onto 
the potential sphere (s) is done by introducing a free electron solution y>i«(e, t r) Yl^r) 
from the potential sphere back to the hard sphere, which joins continuously and differ- 
entiable to the partial wave, 0i?i(e, r R ), at s R and continuously to the screened spherical 
wave at a R i. The radial part of this backwards extrapolated free-electron solution, after 
normalizing it to one at its a-sphere, is given by 

V&fer) = ^ (6,r) = fk(K,r) + g a m (n,r) D a m (e), (15) 

<PRl{t,ClRl) 

where D Rl (e) is the logarithmic derivative of <fiRi{e, r) calculated at the hard sphere am. 
This can be determined from the matching condition between (p RL (e,r) and y9 Ri (e, r) 
at r R = s R , 

D Rl (e) = V{<p%ba a )} = -IM^- Z\t Rl ^ SR \]'lY am ^ SR l] - (16) 
riV) WriV, RUt fafaadVitofasdy-Vifafasn)} 1 ; 

The relation between the values of the free electron function at a and the partial 
wave at s is 

^m{e,a R i) = (fRi(e,a R i) = 1 P{0m( e » s r)} ~ v {9ri( k , s r)} ^ 
</>m{e,s R ) fm{e,s R ) f R i{K,s R )V{f% l (K,s R )}-V{g Rl (K,s R )y 

In Fig. 1 we have plotted the logarithmic derivative at a = 0.7 w, where w denotes 
the average Wigner-Seitz radius, and the normalization function y?ra(e, am) given in (|17D 
in the case of fee Ga. The logarithmic derivative is a never increasing function of energy 
and it has a pole above the top of the respective band. Between these poles D m (e) is 
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smooth functions of energy, which varies more slowly than V{<pRi(e, s)}, because a < s. 
The poles of D Rl (e) depend on the representation (a) and they are not related directly 
to the band structure. The <£>j«(e, Or) from the figure was obtained for partial waves 
normalized in the w-sphere. It is always a smooth function of the energy and vanishes 
at the poles of D Rl (e). 

The slope matrix, Eq. (|13D, the logarithmic derivative, Eq. (|16D, and the normaliza- 
tion function, Eq. flTTD, play a central role in the present implementation of the EMTO 
Theory. 



B. Kink cancellation equation 



Using the free electron solutions from fll^D and (|T5| ) and the partial waves 4>Ri(e, r R ) 
we can introduce a complete basis set defined in the whole space. These exact muffin 
tin orbitals or kinked partial waves may be written in the form 

r RL (e,r R ) = (cP a m^r R ) - <f a m (e,r R )) Y L (f R ) + ^(k 2 ,^), (18) 

where the radial part of the functions <p a Rl and ip Rl are truncated outside the sphere of 
radius sr and outside sr and inside ctR , respectively. Moreover, the I < k ow projection 
of the i/j ri function is truncated inside the sphere of radius cir, while the high-/ com- 
ponents penetrate into the hard spheres. The ip RL (e,rR) functions are continuous and 
different iable in the whole space, except at the hard spheres, where they have non zero 
kinks. In Eq. (|T^) the partial waves are renormalized according to Eq. (|T^|) 

$RL\^ r R) = 7 v ( 19 ) 

From Eq. ( jl7|) and ( |I9| ) it is immediately seen that the multiplicativ normalization of 
the partial waves does not enter in the expression of the kinked partial wave. Forming 
a linear combination of the kinked partial waves, 

*i(r) = E $hfa>*R)v%Lj, (20) 

RL 

and asking for the kinks be canceled we arrive to the kink cancellation or screened KKR 
equations 

E ^R'L'RL^j) v RL,j = 
RL 

E a R> [S%vbl($) - S R'R S vl D RL (e 3 )] v RLj = for all R'L'. (21) 

RL 

Here we have /, /' < li ow . The solutions of this equation are the one-electron energies 
and eigenfunctions, which, using Eq. ([5]) are given by 

URL, = y (22) 

^Ri{ej,a R ) 
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It is worth to note that in the final expression of the wavefunction ^(r), Eq. (|5]), the 
backwards extrapolated free electron solution does not enters. 

In the case of translation symmetry in Eq. ( pT| ) R and R! run over the atoms in 
the primitive cell only, and the slope matrix, and thus the kink matrix K R , L , RL as 
well, depend on the Bloch vector k from the first Brillouin zone. In Fig. 2 we plotted 
the diagonal elements of the fee slope matrix (symbols) calculated at the center of 
the Brillouin zone as a function of the dimensionless energy parameter (nw) 2 . In this 
calculation the matrix inversion in (|13|) was performed in real space for 5 coordination 
shells plus the central site using the s,p and d orbitals and 0.7w for the hard sphere 
radius. The figure demonstrates the weak and smooth energy dependence of the slope 
matrix up to the bottom of the continuum, (kw) 2 ~ 6. Therefore in the practical 
solution of the kink cancellation equation ( pT| ) the slope matrix can be estimated using 
a Taylor expansion around a fixed energy k,q, 

Sr 1 l> rl(. k2 ) = ^r'l'rlI^o) + yy^fl'z/RL^oX^ 2 ~~ K o) + (23) 

where the overdot indicates energy derivative. The first and higher order energy deriva- 
tives are calculated analytically as described in Ref . || . In equation fl23|) k 2 is a complex 
energy not too far from Kq. In Fig. 2 the solid lines were calculated with a fourth order 
expansion around k = 0. As one can observe, this expansion gives highly accurate en- 
ergy dependence of the slope matrix over an energy range of approximately (—1, +l)Ry. 



C. The electron density 

In order to construct the new one-electron potential for a self-consistent calculation 
first we need to construct the electron density given by 

n(r) = eV.WI 2 , (24) 
j 

where the sum runs over the one-electron states below the Fermi level e^. In the present 
implementation of the method instead of calculating explicitly the wave functions (|5|) 
and performing the summation in Eq. ( ^4|) we introduce the path operator g R , L , RL (z , k) 
defined for a complex energy z and Bloch vector k by 

K R'L'R"L"( z , k ) 9r" L "rl{ z ^) = Sr'rSl'l- (25) 

R"L" 

This function is analytic in the complex plane and it has poles at the one-electron 
energies along the real axis. Therefore, using the residue theorem, for the total number 
of electrons we find 

N(e F ) = -L I I g% LIRL {^)K a RLR , L ,{z^)d\,dz, (26) 

Z7U Je F R , VRL JBZ 

where the first integration is performed on a complex contour and the second one in 
the first Brillouin zone. The contour is chosen in a way that it cuts the real axis below 
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the bottom of the valence band and at 6f- In fl26|) 1,1' < k ow . The z dependent partial 
waves, 4>Ri(z,r R ), and logarithmic derivatives, D Rl (z), are obtained by solving Eq. @ 
for complex energy. The energy derivative of the kink matrix, 



Kr'L'rl{ z -> k) 



OR' 



Swl'rl{z- Vq^) - 5rir5vl D a RL {z) 



(27) 



is calculated by taking the derivatives of Eq. ( ]T5| ) and (|2"3"|), where the energy derivatives 
of the basis functions {f a ,g a } are calculated analytically. The energy derivative of the 
logarithmic derivative function is given by M 



dV{<p m (z, s R )} 



Io R <t>Ri(z,r R ) r\ dr R 



dz s R (p Rl (z,s R 
Because the eigenvectors are normalized as (see Ref. 0) 

#*(r) ^-(r) dr -- 



E 

R'L'RL 



V R*L',j R'L'RL (^) 



J RL,j 



(28) 



(29) 



the expression fl26|) gives the exact number of states at the Fermi level for the muffin 
tin potential (g). In ( |2"9"D the negligible terms due to the overlap between s-spheres are 
omitted 0. 

Inside the unit cell at R the electron density in terms of the path operator can be 
expressed as 



n(r R 



2^ I E Z RL'( Z > r n) J BZ 9rl'rl(^ k)dkZ a RL (z, r R )dz, 



(30) 



EF L'L 

where we have introduced the functions 



Z% 



til 



[z, r R ) Y L (r R ) if I < l low and r R < s R 



RL\ z i r RJ 



^Rl( z ^ r R) Yl{tr) 

-ji{z,r R ) Y L (r R ) 



if I < how and r R > s R 
if I > k ow for all r R 



(31) 



and where the sums over /' and / include the high-/ terms as well. These functions are 
equivalent to the scattering solutions of Faulkner and Stocks [|TU|. In Eq. (^) we have 
introduced the following matrix 



9rL'Rl( z ,^-) 



9rl'rl( z ^) ^ ^ I' — how 

J2r"L" 9ru r" L''( z ityS R ii L ii rl (z,\s) if V < li ow and I > li ow 

T1r"L" S R LiRiiL'r(z,]s.)g R ,, L ,, RL (z,]i) if I' > li ow and I < li ow 
J2r"l" J2r"'l"' Srl'R"L"( z i k) 

, x 9r"L"R"'L">(, z i k)S RI „ L „ IRL (z, k) if l,l>li ow 



(32) 



where the high-low and the low-high subblocks of the slope matrix are calculated by 
the usual blowing-up technique |[I . 



In principle Eq. (g6|) and (p0|) give the exact number of states and electron density. 
However, in some cases, like for the metals from the II B and III — V A groups, where 
one of the d bands is completely filled, around the top of this band the normalization 
function ( |17D goes through zero. This happens, for example, in the case of fee Ga 



S 



around the top of the 3c? band, as it can be seen from Fig. 1. For this energy not only 
the logarithmic derivative but also its energy derivative D Rl , appearing in the diagonal 
of the K R , L , RL matrix, has poles. In order to cancel these nonphysical poles we rewrite 
the expression for the number of states as 



N{e F ) 



— j> G(z)dz, 



(33) 



where 



R'L'RL 



G ( z ) = E lg a R /LI R L (z,k)k a RLR/L ,(z,k)dk-Y: 



BZ 



RL 



D a R i( z ) 
D a Rl {z) 



^7 

e D Z 
Rl 



'Rl 



(34) 



with 1,1' < l[ ow , and that of the electron density as 

™M = 77- i Y, Z RL'{ Z ^R) I g a RL>RL( Z , k ) dkZ RL( Z , r R) dz 
2,1X1 Je F L , L JBZ 



2m Je F Y a-R D h( 



EE 



Z RL\ e Rli r R 



(35) 



e Rl 



'Rl\ c RlJ 



where e Rl are the zeros of the logarithmic derivative function, D Rl (e). Because the 



logarithmic derivative is a smooth decreasing function of energy e^'s can be easily 
determined with high accuracy. The second and third terms from the right hand side of 
(j35|) are included only for I < k ow . Using the fact that the residuum of the l/D Rl around 
e Rl is l/D Rl , it is easy to show that the poles of D R i{z) and those of l/ifRi(z,an) 2 are 



canceled out in Eqs. (p4|) and (pq). 

From the electron density (|35|) we determine the overlapping muffin tin wells and 
repeat the iterations until self consistency, of the total energy for example, is achieved. 
In the present implementation of the EMTO Theory the solution of the Poisson's equa- 
tion involves the SCA for the shape of the Wigner-Seitz cell, therefore, the construction 
of the muffin tin potential will be discussed only within this context. 



III. THE SCA-EMTO METHOD 

Equations (|21|) and (|33| - p5D , derived in the previous section, constitute the basis of 
the present method. In order to perform a self-consistent calculation one constructs the 
electron density from the solutions of the kink cancellation equation and calculates the 
new one-electron potential. In this section we describe these steps using the SCA for 
the shape of the Wigner-Seitz cell. 

In the SCA, for solving the Poisson's equation, we substitute the Wigner-Seitz cells 
by spherical cells with volumes equal to the volumes of the real cells. If Qr denotes 
the volume of the Wigner-Seitz cell (Voronoi polyhedron) centered at R we have = 
Q WR = y w ri where wr is the atomic sphere radius. Thus within the SCA, like in the 
conventional ASA, the whole space is "covered" by the Q WR spheres. 
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A. The SCA charge density 



During the self-consistent calculation the Fermi level of a N electron system is 
determined by solving the N(ep) = N equation, where N(ep) is given in fl33|). For this 
€f the electron density is constructed from Eq. (|35|) . Due to the normalization (|29|) the 
so constructed density is exactly normalized within the unit cell cell but not within the 
SCA spheres of volumes £l WR . Therefore, in order to solve the Poisson's equation within 
the SCA we have to renormalize the total density inside the spheres. In the present 
implementation of the method this is realized by 

n SCA (v R ) = n(r R ) + aY 00 (r R ), (36) 

where the site independent a constant is determined from the condition of the charge 
neutrality within the whole unit cell 

W n SCA {v R )dv R = Y,Z R . (37) 

R Jn ™ R R 

Here Z R denotes the nuclear charge at R. The sum runs over the atoms from the unit 
cell, and the integrals are performed inside the SCA spheres. Throughout this section 
the charge density is normalized within the SCA spheres according to ( |3"6"D and (P7|). 
however, for the sake of simplicity we neglect the SCA index for the n SCA {r). 



B. The SCA muffin tin potential 

The spherical symmetric potentials, v R (r R ), that enter in Eq. (|3|) have to be chosen 
in a way that, together with the parameter v , to give the best approximation to the 
full potential v (r). The original idea in Ref. @ is to minimize the mean of the squared 
deviation between the left and the right hand side of Eq. @. This leads to a set of 
integral or differential equations for v R (r R ) and Vq. In the non-overlapping muffin tins 
case the equation for v R {r R ) reduce to the well known expression 

v R (r R ) = — J v(r)dr R , (38) 

and vq reduces to the muffin tin zero, i.e. to the average of the full potential calculated 
in the interstitial region, 

R R 6 

where Q is the volume of the region where the approximation (0) is valid (unit cell), 
and V R denotes the volume of the potential sphere. 

In the overlapping muffin tins case the equation for the vq can be written in the 
following simple form 

Y,J[J K(^) - M r R dr + v o = ^y n v(r)dr, (39) 
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while the equation for v R (r R ) involves terms coming from the overlapping region, and 
which give rise to kinks of Vr{tr) when tr touches other muffin tin spheres. In the 
present implementation of the method, instead of solving the vr(tr) equations, we all 
the time, for non-overlapping and for overlapping muffin tin wells as well, fix the v R (r R ) 
functions to the spherical average of the full potential given by (|38|). In this case from 
Eq. (p9D we get the expression for the vq as 



/ v(r)dr — / v(r)dr 
Jn R Jv R 



or 



\S n i v(r)dr - / n „ v(r)dr 



(40) 



R 



(41) 



where Q R is the real interstitial within a Wigner-Seitz cell centered at R with volume 
Qr, and Qft is that part of the potential sphere that is outside of the cell Q R , i.e. 



Q = Mr and Q R - Q°r" =Q r - V r . 



(42) 



Eq. ( PD assumes the knowledge of the full potential v(r) in Q R and Q R regions. How- 
ever, the time consuming calculation of the full potential can be avoided by using the 
SCA for the unit cell. In the non-overlapping SCA case, i.e. sr < wr, we have 



47T 



and 



v(r)dr 



v(r)dr 



r R dr R , 



v(r)dr R 



r R dr R , 



(43) 



while for the overlapping SCA case, i.e. sr > Wr, we have 

- n£ = -4tt 

and 



r 2 R dr R , 



v(r)dr 



v(r)dr 



"R 



W R 



v(r)dr R 



r R dr R . 



(44) 



From these equations we get the expression for the parameter vq valid within the SCA 



E 

R 



s R 



R 



v(r)dr R 



dr R /J2 Wr 



(45) 



R 



where W R = 47r(w R —s R )/3. Therefore both of the vrItr) function and the Vq parameter 
are given in terms of the spherical symmetric part of the full potential. 
The many-body part, fi xc [n(r)], of the one-electron effective potential, 
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v(r) = v c (r) + ^[n(r)], (46) 

is calculated within the local density or generalized gradient approximation, while the 
electrostatic part is derived solving the Poisson's equation, 



V 2 v c (r) = -8tt 



n(r) - Z R 5(r R )}, (47) 

R 



for the electronic and nuclear charge densities. The electrostatic potential can be di- 
vided into intercell and intercell component. The spherical symmetric part of the in- 
tercell or Madelung potential is given by 

Vr(tr) = - E M rlr , l ,Q r , l , with L = (0,0), (48) 

W R'L' 

where M R l R 'l' is the Madelung matrix, which can be evaluated by the usual Ewald 
technique, and 



Qrl = 7^J~~T I \n R (r R ) - Z R 5(r R )} Y L (f R ) dv R . (49) 

21 + 1 Jvl wr K w ' L J 

The Hartree part of the intracell Coulomb potential can be obtained as the solution 
of the Poisson's equation using the proper boundary condition at the atomic sphere 
radius. Alternatively, this term is given by 



v^Vr) = I 67T \A^ Rr ' R nR ( r '^ dr R + f™R r R n R( r R) dr R\ for r R< w R ( 5Q j 
871--^- Sq R r' R n R {r' R )dr' R for r R >w R ' 



that is valid inside the potential sphere s R , for s R > w R as well as for s R < w R . The 
total potential within the potential sphere is obtained as the sum of Eq. (pE8D, (|50|). 
the Coulomb potential of the nucleus and the spherical symmetric exchange-correlation 
potential, namely 

2Z 

vr{t r ) = v R + v R (r R ) h ^ xcR (r R ). (51) 

rR 

If the spherical symmetric part of the exchange-correlation potential is approximated 
by HxcRinRirR,)] besides the higher order multipole moments from (f4"8"D, which in many 
cases can be neglected, all of the potential components from ( |5T| ) depend only on the 
spherical symmetric density n R (r R ). 

Within the SCA-EMTO method the atomic and potential spheres can be and usually 
they are chosen differently. The sizes of the atomic spheres, w R , are fixed by the volume, 
and the ratio between them should be chosen in a way that minimizes the errors coming 
from approximate solution of the Poisson's equation. We have found that the best 
representation of the potential can be achieved by choosing the potential sphere radii, 
s R , larger or equal with the atomic sphere radii. For an optimal choice of the potential 
spheres the potentials at s R should be the same, i.e. v R (s R ) ~ constant for each R, 
and this constant should have the maximum possible value for linear overlaps bellow 
30 - 40%. 
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IV. APPLICATIONS: TEST CALCULATIONS 



In this section we present a few applications of the SCA-EMTO method. We chose 
particular systems where the conventional ASA based methods failed and the inclusion 
of the correction terms or of the exact potential seemed to be unavoidable. First we 
describe the most important numerical details and after we analyze the present results 
comparing to the available full-potential calculations. 



A. Numerical details 

The hard sphere radii are chosen at a R = 0.7w. In the matrix inversion from Eq. 
([□D we includ 79 sites in the case of fee-based structures (fee, Ll 2 and Llo), and 89 
sites in the case of free-based structures (bcc,B2 and B32). The Taylor expansion of 
the slope matrix is carried out for k = and includes terms up to the fourth to sixth 
order energy derivative. 

The path operator is calculated for 16 — 32 complex energy points, depending on 
the band structure, distributed exponentially on a semi-circular contour. The fc-point 
sampling is performed on a uniform grid in the 3D Brillouin zones. All the calculations 
are scalar-relativistic and employ the frozen-core approximation. The basis set include 
s, p and d orbitals and the valence electrons are treated self-consistently within the 
local density approximation to density functional theory using the Ceperley and Alder 
T2|| exchange-correlation functional and parametrized by Perdew and Wang fl3 |. The 



atomic sphere radii, w R -s, are chosen in a way that the atomic spheres should have the 
same volumes as the corresponding Voronoi polyhedra. The electrostatic and exchange- 
correlation contribution to the total energy is calculated within the SCA as described, 
for example, in Refs. [|Il] , |i"5| . The kinetic energy is given by 

] p p w R 

— & zG(z)dz - V / v R (r R )n R (r R )r R dr R , (52) 

llTt Je F p Jo 



iSCA _ 

2%i 



where the first term from the right hand side is the sum of the one electron energies 
and G(z) is given in ([34]) . 



B. Results 

Before starting on the evaluation of the results we address the question of the accu- 
racy of the Taylor expansion for the slope matrix, Eq. ([£J). In Fig. 3 the total energy 
of fee Cu is shown for different number of terms included in the Taylor expansion. 
The inclusion of the fourth order energy derivative term changes the total energy by 
1.13 mRy. The effect of the fifth order term is already less than 0.2 mRy and that of 
the sixth order term is about 0.04 mRy. Therefore, we conclude that for a reasonable 
accuracy it is sufficient to include five terms in the expansion of the slope matrix, viz. 
up to the fourth order energy derivatives. In the case of open structures, wide energy 
bands or semicore states, however, more terms should be included . 
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The variations of the SCA-EMTO total energy with the potential sphere radius for 
fee and bee Cu are shown on Fig. 4. For this test calculation the total potential from 
( fj5|) was weighted by the valence part of the total density according to 

Efi Cj [In(v)v(r)dv R } dr R 
V ° ' Er ir R R r R [Jn(v)dr R ] dr R ' ^ 



and the cut-off in (p5|) for the spherical part of the total density was l max = l' max = 8. 
The inscribed sphere radii are s^ cc = 0.91w/ cc and s l bcc = 0.88wb cc , where the theoretical 
atomic sphere radii, Wf cc and Wb cc are shown on the figure. The total energy in both 
of the fee and bee structures beginning from s ~ 0.80w becomes almost flat with a 
negligible slope up to s ~ 1.20w, which means about 32 — 36% linear overlap for the 
fee and bee structures, respectively. For s > 0.80k> the further increase of the potential 
sphere radius has little effect on the energy, that means the potential in the corners of 
the Wigner-Seitz cell is, with a very good approximation, constant. However, for big 
overlaps, s > 1.20k;, the errors coming from the overlap region, and neglected in the 
kink-cancellation equation and in the charge density as well, become important ||. 

There is a comprehensive study of the structural stability of the transition metals 
done either by full-potential or by muffin-tin or ASA based methods. In the latter 
case correction terms are needed |17| for calculation of the accurate total energies. The 



conventional ASA without correction terms gives, for example, with about 2 mRy lower 
total energy for the Cu in the bee phase than in the fee phase. This underestimation of 
the bee total energy is due to the incorrect kinetic energy term, and the inclusion of the 
exact Hartree energy would lower even more the bee energy. From a more sophisticated 



full-potential method |TB[ a structural energy difference of Eb cc — Ef cc « 0.5 mRy was 
obtained. This number should be compared with our results of 0.4 mRy from Fig. 
4. One should note that this difference is almost constant for a wide range of linear 
overlap. 

The second example is the binary \ji x A\i_ x ordered compound in different phases. 
There are two reasons of this choice: i) this system is very well studied through ac- 



curate full-potential calculations |19| , |20[| , and ii) most of the experimentally observed 
interesting trends, the contraction of the volume of the Al-based alloys, the asymmetric 
heats of formation with respect to the equiatomic concentration etc., can not be re- 
produced by an ordinary spherically symmetric calculation. Besides the three different 
compositions, x = 0.25, 0.50, 0.75, we consider the pure Al and Li limits in fee and bee 
phases as well. For x = 0.25 and 0.75 the calculations were performed only for the LI2 
structure, while for x = 0.50 we considered three different structures: LIq, B2 and 532. 

In Fig. 5 the charge density contour plots are shown for pure fee Al, and Al 3 Li in 
Ll structure, as calculated from Eq. (^) using s,p,d basis set, and the maximum 
orbital quantum number /' included in (|T2"D was 10. The agreement between these plots 
and those from Ref. |^T| is very good. 



The calculated equilibrium Wigner-Seitz radii and the bulk moduli are tabulated 
in Table I and plotted in Figs. 6 and 7. In the case of pure Al and Li the structure 
energy differences and in the case of compounds the heats of formation are included in 
the Table and plotted in Fig. 8. The heat of formation is defined as 

AH = E LixAh _ x -xE u -(l- x) E Ah (54) 
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where all the energies are obtained for the proper equilibrium volume and they are 
expressed per atom. In Figs. 5-7 and Table I the full-potential values from Ref. [20] are 
also included. 

The mean deviations between the present and the full-potential results for the equi- 
librium radii, bulk moduli and heats of formations are 9.8%, 7.5% and 17%, respectively. 
Taking into account the minor discrepancies between the numerical details used in the 
calculations, for example the exchange-correlation functional, the way the core elec- 
trons were treated etc., and the fact that the full-potential methods have their own 
error limits as well, we can conclude that the agreement between the two sets of results 
is very good. One should appreciate how well the trends obtained in the full-potential 
calculation are reproduced by the present method. 



V. CONCLUSIONS 

We have presented a self-consistent implementation based on the Green's function 
technique of the Exact Muffin-Tin Orbitals Theory, developed by O.K. Andersen et 
al. The accuracy of the present implementation was tested on different systems, where 
we have found a good agreement between the present results and the results obtained 
by full-potential techniques. In order to gain some experience about the efficiency of 
the present method we compare the CPU times of a self-consistent calculation of the 



tight-binding ASA-LMTO method |TT|, based also on the Green's function technique, 
and that of the SCA-EMTO method. We found that the present implementation of the 
SCA-EMTO method needs with about 3 times larger CPU time than the tight-binding 
ASA-LMTO method. 

Finally we remark that if the radii of the potential spheres are chosen to be equal 
with the radii of the atomic spheres, i.e. w r = sr, the SCA-EMTO method can be 
considered as an ASA based Green's function technique that involves the so called 
combined correction term ||22|| . It gives exact one electron energies and charge densities 



for the optimized overlapping muffin tin wells. The natural extension of the present 
SCA-EMTO method to compute the total energies from the output total charge density 



via the Full Charge Density technique p3,24] is in progress 
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VI. APPENDIX 



During the self-consistent procedure the Eq. (ESI), (P5|), (|35J) and (pTJ) are solved 
iteratively. In order to construct the electron density ([55|), as the input for the next 
iteration, we have to invert the kink matrix for each complex energy z along the contour 
and for each Bloch vector k from the Brillouin zone. For a reasonably high accuracy 
we need at least a few hundreds of Bloch vectors in the irreducible part of the Brillouin 
zone, therefore this solution means the most time-consuming step of the self-consistent 



procedure. Here we apply a similar "two-step" scheme introduced in Ref. [^5j in order 
to reduce the number of time-consuming iterations. Within this scheme after each 
iteration an approximate charge self-consistency is achieved by solving self-consistently 
the following equation written for the k-integrated path operator 



g\z) = [l + g a0 (z) (D a0 (z) - D^z))}' 1 g a0 (z), 



(55) 



where where the index denotes quantities obtained from the previous iteration, and 
which are kept fixed during the solution of Eq. (155]). 

In the expression for the number of state (|33|) , we need the k-integrated trace of 
the product between the path operator and the energy derivative of the kink matrix, 
therefore, a similar equation to fl5"5]) has to be established for the quantity 



Gr>l>Rl( Z ) — 9R'L'RL( Z ik) Kft LR , L ,(z,\i) <lk. 

J BZ 



(56) 



Using the definition of the kink matrix after some manipulations we arrive to the equa- 
tion 



Gr>l'rl{ z ) 



9r'l'rl\ z ) 



Givl>rl( z . 
9r'l'rl( z ) 



+ Wz/l a R (D R ?(z) - D Rl (z] 



(57) 



It is worth to note that in this expression we do not have matrix multiplication. Finally 
we mention that as soon as the self-consistency is achieved, i.e. 



D a m °(z) and D a m (z) 



D a m °(z), 



(5f 



both equations, ( |55j ) and Q57D, become exact. 
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TABLES 



TABLE I. The theoretical atomic radii, bulk moduli, heats of formation and fcc-bcc struc- 
ture energy differences of the ordered AlLi compounds, and pure Al and Li. The present 
SCA-EMTO results are compared to the full-potential FLAPW results, and to the experi- 
mental data. 





Structure 




Present calculation 


Full-potential 1 


Experimental 


Al 3 Li 


Ll 2 


S(Bohr) 


2.904 


2.935 


2.934 3 






B(GPa) 


74.27 


70.31 


66.0 4 






AF(mRy) 


-8.21 


-8.3 




AlLi 


Ll 


S(Bohr) 


2.881 


2.917 








B(GPa) 


51.04 


50.41 








Aif(mRy) 


-11.53 


-10.25 




AlLi 


B2 


S(Bohr) 


2.837 


2.876 








B(GPa) 


55.86 


42.09 








Aif(mRy) 


-13.94 


-10.45 




AlLi 


B32 


S(Bohr) 


2.863 


2.910 


2.928 5 






B(GPa) 


59.91 


57.75 








AH(mRy) 


-21.69 


-16.60 


-18.5 6 


AlLi 3 


Ll 2 


S(Bohr) 


2.932 


2.965 








B(GPa) 


31.45 


28.37 








A#(mRy) 


-6.71 


-5.0 




Al 


fee 


S(Bohr) 


2.920 


2.946 


2.991 2 






B(GPa) 


91.76 


82.20 


72. 8 2 






E bcc - Ef cc (mRy) 


2.90 


4.59 




Al 


bec 


S(Bohr) 


2.930 


2.951 








B(GPa) 


84.46 


84.18 




Li 


fee 


S(Bohr) 


3.112 


3.124 








B(GPa) 


15.47 


13.64 








E bcc - Ef cc (mRy) 


0.29 


0.50 




Li 


bec 


S(Bohr) 


3.117 


3.128 


3.237 2 






B(GPa) 


15.64 


15.25 


12. 6 2 



1 FLAPW calculation, Ref. 2 Experimental, Ref. p6 
3 Experimental, Ref. [27]; 4 Experimental, Ref. 



Experimental, Ref. | 29| ; Experimental, Ref. [30] 
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FIGURES 



FIG. 1. The energy dependence of the logarithmic derivative D^(e) (solid line), and 
normalization function l -Pm{ e i a R) (dashed line) for the fee Ga. 

FIG. 2. The diagonal elements of the fee slope matrix in the k = (0,0,0) point from the 
Brillouin zone versus (k w) 2 . The numbers in parenthesis denote the (I, m) quantum numbers. 
The Taylor expansion included terms up to the 4th order energy derivative. 

FIG. 3. Self-consistent SCA-EMTO total energy of fee Cu for different number of terms 



included in the Taylor expansion for the slope matrix, Eq. (|23|). 



FIG. 4. Total energy versus potential sphere radius, s, for the fee and bee Cu. The 
calculations were done at the theoretical equilibrium atomic sphere radii shown in the figure. 

FIG. 5. Charge density contour plots for ALjLi in LI2 structure and fee Al in units of 
0.01 electrons/ Bohr 3 obtained as the output of a self-consistent SCA-EMTO calculations. 

FIG. 6. The change of the atomic radii of the ordered AlLi compounds relative to the 
radius of the fee Al. The open symbols show the results obtained by the full-potential FLAPW 
calculation from Ref. ^(|. The present SCA-EMTO results are shown by closed symbols. The 
lines connect the results for the /cc-based structures. 

FIG. 7. The theoretical bulk moduli for the ordered AlLi compounds. For the notation 
see caption of Fig. 6. 

FIG. 8. The theoretical heat of formations for the ordered AlLi compounds. For the pure 
Al and Li the fcc-bcc structure energy difference is shown. For the notation see caption of Fig. 
6. 
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Energy (Ry) 




Energy *w 2 




Number of terms in the expansion 



Linear overlap, fcc/bcc (% 



-0.82 



2/4 



18/22 



CD 
C 
LU 



-0.83 - 



-0.84 - 



-0.85 - 



fee (w = 2.613 Bohr) 



bee (w = 2.614 Bohr) 




35/39 



/ _ 




-0.86 - 



1.6 



2.0 



2.4 



2.8 



3.2 



3.6 



Potential sphere radius (Bohr) 



(a) AI3Li (100) 




(b) AI3Li (111) 




(C) Al (111) 
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